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ABSTRACT 


This report discusses the Integration of the vector differential equation 
x * f(x, t) from time tj to t| + j, where only the values of are available for 
the integration. No previous values of x or x are used. Using an orbit inte- 
gration problem, comparisons are made between Taylor series integrators and 
various types and orders of Runge-Kutta Integrators. A very outstanding fourth 
order Runge-Kutta type integrator for orbital work Is presented. Approximate 
(there may be no exact) fifth order Runge-Kutta integrators are discussed. Also 
discussed and compared Is a self starting Integrator using af/ax. A numerical 
method for controlling the accuracy of integration Is given. And, the special 
equations for accurately Integrating accelerometer data are shown. 



1 . INTRODUCTION 


Within recent years, the increasing use of sequential data processors 
(e.g. Kalman filters) has increased the need for high quality, self-starting 
integrators. A self-starting integrator is one which takes a current estimate 
of the state vector, x^, and, using the equation x * f(x^ t), propagates this 
estimate ahead aT seconds to obtain x^. Non self-starting integrators (e.g. 
Adams-Moulten integrators) need the previous values of x^, x. J, x^.g* etc * in 
order to obtain an estimate of x^ +1 . 

This report will discuss several self-starting integrators. An empirical 
evaluation of these integrators will then be made with regard to their efficacy 
in an orbit determination problem. A numerical method of determining the ac- 
curacy of these integrators over each integration step is shown. This determina- 
tion is quite useful in such areas as calculating interplanetary trajectories, 
where equal accuracy integration step sizes may vary from a few seconds (close 
to a planet) to as long as several hours. 

A section has been included which shows how to accurately integrate ac- 
celerometer data. This data generally comes into the integrator as the integral 
of sensed acceleration, not *ensed acceleration itself. 

To my knowledge, no one has ever developed a fifth order Runge-Kutta in- 
tegrator for the vector differential equation ;x * £(x, t). We will show that 
6 evaluations of _f(x, t) are needed, not 5 as one might expect. Further, even 
with 6 evaluations there may be no exact set of integrator constants. Two ap- 
proximate sets of fifth order integrator constants are shown and evaluated. 
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2. INTEGRATION BY TAYLOR SERIES 


One of the oldest and best known methods of numerically Integrating is the 
Taylor series expansion. It is the standard o'. comparison against which all 
other methods are evaluated. Despite the fact that the method is old, doesn't 
mean that it is obsolete. Indeed it is frequently and efficiently used in digital 
computer programs. Integration of the vector differential equation, x * f(x, t), 
is accomplished by 


—1+1 * £j * ij AT + *1 TT + TT* '*i TT + 


( 1 ) 


where 

x = f (x , t) 

3X . 3X 3f 3f 

X = X + = f + 

- 3X - 3t 3X - 3t 


3JC . 3X 
3X - + 3t 


3X 


3X 


X = X + 

- 3X — 3t 


etc. 


The partial derivative of a v-ctor, £, with respect to a vector, z, is defined 
to be the matrix 
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where the subscripts indicate the particular element of the vector in question. 
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3. SECOND ORDER RUNGE-KUTTA TYPE INTEGRATOR 

A second order Integrator is one that gives a solution which agrees with 

2 

a Taylor series expansion up through the aT 6, terns. Such an Integrator is 


■ ATf(Xj, t| + s^aT) 

(2) 

k 2 * ATf(xj + a^k}* t^ + 6 2 aT) 

(3) 

x.i+1 s + b-jkj + b 2 k 2 + 0(AT 3 ) 

(4) 

The integration constants in the above equation are not arbitrary. In order 

3 

to achieve an error of order aT , the following constraint equations must be 
satisfied 

b 2 a 1 = 1/2 

(5) 

b 1 + b 2 = 1 

(6) 

b-j 5 1 + b 2 d 2 = 1/2 

(7) 


Note that there are only three equations in five unknowns. The constraint 
equations could be rewritten as 
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6-j and 6 2 can be experimentally chosen to reduce the 0 (aT ) truncation error. 

A special case arises when 6-j s fi 2 . constraint equations for this case are 
shown In Table 2. Here, b 2 Is chosen to reduce the truncation error. 


5 1 *7 


h * 7 


a l 




TABLE 2: 2ND ORDER CONSTRAINTS WHEN 5 1 - fig 


It Is 
be achieved 


interesting to note that fourth order Integration of x 
by adding the two additional constraint equations 


f(t) can 


b l 6 l + b 2 S 2 = 1/3 

b^ 3 + b 2 s 3 = 1/4 (9) 


There are two solutions of the five constraint equations (equations (5) tnrough 
(9)). These solutions are shown in Tables 3 and 4. 

6 1 = (3 -V3)/6 = .2113 2486 5405 1888 

« 2 * (3 + V3)/6 * .7886 7513 4594 8112 

a-j » 1 bj * 1/2 b 2 = 1/2 


TABLE 3: 2/4 ORDER INTEGRATOR CONSTANTS, SET #1 





6, » (3 +^?)/6 » .7886 7513 4594 8112 

« 2 ■ ( 3 -&/* * .2113 2486 5405 1888 

a 1 = 1 b 1 * 1/2 i> 2 = 1/2 


TABLE 4: 2/4 ORDER INTEGRATOR CONSTANTS, SET #2 


When x = f(t), the constants in Tables 3 and 4 will cause the solution 

to be 


-i+1 


means approximate value of x^ + -j • A similar equation also applies to any 
element of fjx, t) which is a function of time alone. Note that the general 
solution of x = f(x, t) is still only second order, hence the terminology 2/4 
order. 

If, instead of using equation (9) as an additional constraint equation, 
we use the equation 


6 -| — 0 ( 10 ) 

then we obtain a third order solution when x = f(t). The 2/3 order integrator 
constants are shown in Table 5. 


fi 1 « 0 6 2 = 2/3 a ] = 2/3 b ] » 1/4 b 2 * 3/4 


TABLE 5: 2/3 ORDER INTEGRATOR CONSTANTS 
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4. THIRD ORDER RUNGE-KUTTA TYPE INTEGRATOR 


Third order Integration of x ■ f(x, t) can be achieved by 


jtj * AT ffx^, t^ + <5.|AT) 

(11) 

kg * aT f(x. + a-jkp t| + 6gAT) 

(12) 

— aT f (x^ ^ b-| k i + bpkg t t^ ^ A^aT) 

03) 

x. 1+ -j = x ^ + c-jJc-j + c 2 Ji2 + + 0(aT 4 ) 

(14) 


The constraint equations for the integration constants are 


o> 

II 

O 

(15) 

a l = 6 2 

(16) 

b l + b 2 = 6 3 

07) 

b 2 C 3 6 2 = 1 / 6 

08) 

C 1 + c 2 + c 3 - 1 

09) 

c 2 6 2 + c 3 6 3 * 

(20) 

c 2 6 2 + c 3 6 3 * 

(21) 
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Here again there are two more unknowns than there are equations. Thus all the 
Integrator constants can be expressed as functions of 6 2 and 63 , as shown In 
Table 6 .* & 2 an< * 6 3 can experimentally adjusted to reduce the 0 (aT 4 ) trunca- 
tion error In equation (14). 



TABLE 6 : THIRD ORDER CONSTRAINT EQUATIONS 


Two special cases arise when <s 2 = <$3 an< ^ when <$3 = 0. The constraint 
equations for these cases are shown in Tables 7 and 8 . In both cases c 3 is 
emperically determined to reduce truncation error. 



TABLE 7: THIRD ORDER CONSTRAINTS WHEN & 2 = 63 


* Note from equation (18) that b 2> C 3 , s 2 * 0 are not allowable solutions. 
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w° 

2 2 
5 2 a 3 Vi 

b ' 

h - 1 

b 2-*cJ 

C 1 3 l * C 3 

c 2 = T 


TABLE 8: THIRD ORDER CONSTRAINTS WHEN fi 3 = 0 


Fifth order integration of x * f(t) can be achieved by adding the two 
additional constraint equations 


2 6 2 + C 3 6 3 = 1/4 

(22) 

2 6 2 + c 3 6 3 = ly ^ 5 

(23) 


There are two solutions of the nine constraint equations (equations (15) 
through (23)). These solutions are shown in Tables 9 and 10. 


6-j = 0 




6 2 = (6 -V6)/10 = .3550 

5102 

5721 

6822 

63 = (6 '+V6)/10 = .8449 

4897 

4278 

3178 

a ] = (6 -yf6)/10 = .3550 

5102 

5721 

6822 

b-j = -(54 + 19^6)/250 * - 

.4021 

6122 0451 5215 

b 2 = 2(51 + 11^/125 = 

.247 

1101 

9472 98393 

c 1 * 1/9 = .1111 1111 1111 1111 


c 2 = (16 +V6’)/36 * .5124 

8582 

6188 

4216 

c 3 = (16 - V?)/36 » .3764 

0306 

2700 

4673 


TABLE 9: 3/5 ORDER INTEGRATOR CONSTANTS, SET NO. 1 
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6 1 * 0 

6 2 » (6 +>/?)/ 10 « .8449 4897 4278 3178 

6 3 - (6 - V ?)/10 = .3550 5102 5721 6822 

a 1 * (6 +V6)/10 « .8449 4897 4278 3178 

b ] * -(54 - 19^6)/25Q * -.02983 8779 5484 7846 

b 2 = 2(51 - llV6)/125 = .3848 8980 5270 1607 

c 1 * 1/9 * .mi mi im nn 

Cp * (16 -yj6)/36 = .3764 0306 2700 4673 
c 3 = (16 +y[6)/36 = .5124 8582 6183 4216 

TABLE 10: 3/5 ORDER INTEGRATOR CONSTANTS, SET NO. 2 


When x = f(t), the constants in the two tables above give a solution of 


x_. jn a X, + X.aT + X^ + ” 


ii+1 - *1 T 


d 5 x. t 5 d 6 x. aT 6 

.99 - 4 - Hr + 

dt b b- dt b 
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5. FOURTH ORDER RU.1GE-KUTTA TYPE INTEGRATOR 


Fourth order integration of x = f(x, t) can be achieved by 

jc-j = AT f (x^ , t. + a^T) 

k 2 = aT f(x. + a-jkp t i + 5 2 aT) 

k_ 3 = aT f(x. + b-jk^ + b 2 Jc 2 , t^ + 5 3 aT) 

JI4 ” + *"1—1 + c 2— 2 + c 3^3 1 ^-j + 4 ^aT) 

-i+1 = *i + d lil + d 2— 2 + d ^3 + d 4^4 + 0( * t5) 
The constraint equations for the integration constants are 

<1=0 

a l = 6 2 

- b l + b 2 = 5 3 

C 1 + c 2 + c 3 = 5 4 

b 2 d 3 5 2 + ^ c 2 5 2 + c 3 d 3^ d 4 ~ 1/,fi 

b 2 d 3 6 2 6 3 + ^ C 2 6 2 + c 3 6 3^ d 4 5 4 = 
b 2 d 3 5 2 + {c 2 d 2 + c 3 5 3 )d 4 = 1/12 

°2 C 3 d 4^2 — 1/24 


(24) 

(25) 

(26) 

(27) 

(28) 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 
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dj + dg + d^ + d^ * 1 

(37) 

d 2 6 2 + d 3 { 3 + d 4 6 4 * 1/2 

(38) 

d 2 «2 + d 3 6 3 + d 4 6 4 " 1/3 

(39) 

Vz * d 3 «f ♦ d 4 a| - 1/4 

(40) 


It Is interesting to note that equations (34) and (35) cannot be derived from 
consideration of the scalar equation x = f(x), as can most of the other con- 
straint equations. Consideration of only the scalar equation leads to a con- 
straint equation which is a linear combination of equations (34) and (35). 

It can be shown that the only value of 6^ which satisfies the constraint 
equations is 6 & * 1. As before, the constraint equations can be rewritten as 
shown in Table 11. 6 2 and can be emperically determined to reduce the 0 (aT ) 
truncation error in equation (28). 



TABLE 11: FOURTH ORDER CONSTRAINT EQUATIONS 
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There Is an interesting array of nonacceptable solutions of the constraint 
equations. Some of these are 


« 2 t 0 

H* i 

~ 4«2 * 46 3 + 662^3 + 3^0 unless 6 2 * 1 
Z &2 - 1^0 unless 63 * 63 or unless 63 * 0 
b 2 |i 0 
c 3 #0 

V° 

d 4 l* 0 

These unacceptable solutions give rise to the three special forms of the 
constraint equations shown in Tables 12, 13, and 14. If d 3 is set to 1/3 in 
Table 12, it wi 1 ! be seen that the standard, fourth-order, Runge-Kutta integra- 
tion constants will be obtained. If d 3 is set equal to (1 +l/^2)/3 the 
Runge-Kutta-Gill integration constants will be obtained. 



TABLE 12: 4th ORDER CONSTRAINT EQUATIONS FOR 6 £ * 63 
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* 0 

«2 “ 1 

*3 * J S 4 * 1 

a^ * 1 

b 1 * 3/8 

b 2 = 1/8 

c, • t-c 2 -c 3 c 2 . 

1 


d l 

d 2 = l • d 4 

d 3 “ 1 

TABLE 13: 

4th ORDER CONSTRAINT EQUATIONS FOR S £ = 1 



TABLE 14: 4th ORDER CONSTRAINT EQUATIONS FOR 63 = 0 


Sixth order integration of x * f(t) can be accomplished by adding the 
two additional constraint equations 


d 2 6 2 + d 3 fi 3 + d 4 * 1/5 

(41) 

d 2^2 + d 3^3 + d 4 = 1/6 

(42) 


There are two solutions of the 14 constraint equations (equations (29) through 
(42)). These solutions are shown in Tables 15 and 16. 
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« 2 * (5 -y[s)nO * .2763 9320 2250 0210 
« 3 - (5 +yfi )/ 10 - .7236 0679 7749 9790 


a 1 * (5 -^5^/10 * .2763 9320 2250 0210 
b 1 = -(5 + 3^5)/20 = -.5854 1019 6624 9685 
b 2 = (3+>/5)/4 = 1.309 0169 J437 49475 
c 1 = -(1 - 5^/?)/4 = 2.545 0849 7187 4737 
c 2 * -(5 + 3^/4 = -2.927 0509 8312 4842 
c 3 = (5 -^5)/2 = 1.381 9660 1125 0105 
d 1 = 1/12 = .08333 3333 3333 3333 
d 2 = 5/12 = .4166 6666 6666 66667 
d 3 = 5/12 = .4166 6666 6666 66667 
d 4 = 1/12 = .08333 3333 3333 3333 

TABLE 15: 4/6 ORDER INTEGRATOR CONSTANTS, SET NO. 1 


5-5 




6j * 0 

«2 * (5 +yJs)no - .7236 0679 7749 9790 
« 3 s (5 « .2763 9320 2250 0210 

64 « 1 

^ * (5 +y[s)/}0 = .7236 0679 7749 9790 
b-, = -(5 - 3^5)/20 = .08541 0196 6249 6845 
b 2 = (3 -y[s)/A = .1909 8300 5625 05256 
c ] ~ -(1 + S^5)/A = -3.045 0849 7187 47373 
c 2 = -(5 - 3^5)/4 = .4270 5098 3124 8423 
c 3 = (5 +Vs')/2 = 3.618 0339 8874 98950 
dj « 1/12 « .08333 3333 3333 3333 


d 2 = 5/12 = .4166 6666 6666 66667 
d 3 * 5/12 = .4166 6666 6666 66667 
d 4 = 1/12 * .08333 3333 3333 3333 

TABLE 16: 4/6 ORDER INTEGRATOR CONSTANTS, SET NO. 2 
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When x » f(t), the constants in the two above tables give a so'Lclon of 



*i + *i 




+ 


d6 -i AT 6 301 
dt® ^ 



+ 


A 4/5 order integrator can be obtained from Table 11 if we set 


« 3 * | (3 - 56 2 )/(l - 26 2 ) 

& 2 is then adjusted to give optimum performance. There are, however, 
nonallowable values of $ 2 . They are 

<$ 2 / 0 , . 4 , . 5 , . 6 , 1 , (6 +\ fi ) 


It is sometimes of value to have an estimate of x at t. + 6 aT. 
equations shown in Appendix A, it is easily seen that 


-i+6 = i f + e 1 k 1 * S 2 k 2 * Bjkj + 8^ + 0(iT 4 ) 

where 

s 2 (c 2 6 2 + c 3 6 3 )(2<5 - 36 2 ) - (1 - <5 2 )fi 
3 3 = ~6 6 3 (6 3 - 6 2 ^ c 2 6 2 + c 3 6 3 } ‘ b 2 6 2^ “ 

fi 2 6 3 (6 3 - & 2 ) S - b 2 & 2 (2& - 3S 2 ) 

6 4 * V 6 3 U 3 - & 2 )Lc 2 & 2 + c 3 d 3 ) - b 2 6 2 0 - & 2 ) 


several 


From the 


(43) 


(44) 


(45) 
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(46) 


*2 * ( J ‘ 6 3 4 3 * e 4> 

8^*6“62”®3“ 6^ (47) 

The k's are the same as those used to obtain x^. Thus no new derivative 
evaluations are necessary- 


5-8 



6. FIFTH ORDER RUHGE-KLJTTA TYPE INTEGRATOR 
WITH FIVE DERIVATIVE EVALUATIONS 


Fifth order integration of x = f(£, t) can be attempted by 


-1 

■f(x^» t i + 5 -jAT) 

(48) 

—2 = 

ATf(x. "*■ 3 1 , t^ + 5 2 aT) 

(49) 

-3 

' Tf(x^ *■ b^kj + b 2-2* t i + 5 3‘-"^ 

(50) 

^4 = 

■ Tf (x_. + c 1 k 1 + c 2 k 2 + c-jk-j, t. + i^AT) 

(51) 

II 

-T 

"T£.(Jf. + d-j k_-j + ^ 2—2 + d 3^3 + d 4^4’ *i + 6 5 &T ) 

(52) 

-i+1 

= + e 1 k 1 + e 2 k 2 + + e^ + 0 (aT ? ) 

(53) 

constraint equations for the integration constants are 


5 1 

0 

(54) 

a l * 

'2 

(55) 

b l + 

b 2 ' S 3 

(56) 

C 1 + 

c 2 + c 3 ' '4 

(57) 

d l + 

W V =5 

(58) 

e l + 

e 2 - e 3 + e 4 + e S = 1 

(59) 

e 2 2 

+ e 3 d 3 + + e 5 -: 5 = 1/2 

(60) 
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e 2 6 2 + e 3 6 3 + e 4 6 4 + e 5 6 5 “ 1/3 (61) 

e 2 6 2 + e 3 6 3 + e 4 d 4 + ®5 d 5 “ V 4 (62) 

e 2 6 2 + e 3 { 3 + e 4 6 4 + e 5 d 5 “ 1/5 (63) 

b 2 d 2 e 3 ^ ^2^2 ^ c 3 d 3^ e 4 ^ C^£ d 2 ^3 d 3 ^ ^4 d 4^®5 8 (64) 

b 2 6 2 e 3 6 3 + (c 2 6 2 + c 3 6 3 )e 4 5 4 + (d 2 d 2 + d 3 6 3 + d 4 6 4 )e 5 6 5 = 1/8 ( 65 ) 

b 2 d 2 e 3 6 3 + ^ c 2 6 2 + c 3 6 3* e 4 6 4 + ( d 2 6 2 + d 3 6 3 + d 4 6 4^ e 5 6 5 = 1/10 (66) 

( b 2 6 2 ) 2e 3 + (c 2 6 2 + c 3 6 3 ) 2 e 4 + (d 2 6 2 + d 3 fi 3 + <* 4 « 4 ) 2 e 5 * 1/20 (67) 

b 2 6 2 e 3 + (c 2 6 2 + c 3 d 3 )e 4 + (d 2 6 2 + d 3 6 3 + d 4 6 4 )e 5 * 1/12 (68) 

b 2 d 2 e 3 d 3 + {c 2 6 2 + c 3 6 3 )e 4 6 4 + (d 2 6 2 + d 3 6 3 + d 4 6 4 )e 5 6 5 * 1/15 (69) 

b 2 £ 2 e 3 + (c 2 6 2 + c 3 6 3 )e 4 + (d 2 d 2 + d 3 6 3 + d 4 d 4 )e 5 * 1/20 (70) 

b 2 d 2 c 3 e 4 d 4 + ^ b 2 d 2 d 3 + ^2 d 2 + c 3 d 3^4^®5 d 5 * V30 (71) 

b 2 d 2 c 3 e 4 + tb 2 6 2 d 3 + {c 2 6 2 + c 3 d 3 )d 4 Je 5 * 1/24 (72) 

b 2 6 2 c 3 e 4 6 3 + tb 2 6 2 d 3 6 3 + (c 2 d 2 + c 3 6 3 )d 4 6 4 Je 5 = 1/40 (73) 

b 2 6 2 c 3 e 4 + t b 2 5 2 d 3 + (c 2 d 2 + c 3 6 3 )d 4 ]e 5 * 1/60 (74) 

' VI 20 (75) 
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Since the above equations are so difficult to derive, we mention the 
fact that they have been triple checked for accuracy and may be used with 
confidence by the reader. There are 22 equations In 20 unknowns, leading one 
to suspect that there may be no exact solution of these equations. However, 
the equations are nonlinear and there may be a possibility of a solution. In 
Appendix B we explore this possibility and arrive at the conclusion that there 
Is no exact solution of these equations. Thus it appears that a fifth order 
Runge-Kutta Integrator will require at least 6 derivative evaluations, not 5 
as one would expect. 

We have made though an approximate solution of the above equations. 
Equations (54) through (59) are solved exactly. Equations (60) through (75) 
were solved to an average lo error of + 1.6*1 0“ 6 . The constants are shown in 
Table 17. 


5 1 = 

0 

d 2 = .0013290 82957 

5 3 = .37046 68626 

6 4 

.74985 56689 

<$ 5 = 1.0000 07268 

a l = 

.0013290 82957 



b, = 

i 

-51.257 65184 

b 2 = 51.628 11870 

C 1 = 

165.34 07125 

c 2 = -165.94 52689 

r 3 = 1.3544 12081 

V 

-598.90 49520 

d 2 = 601.87 06163 

d 3 * -3.0940 32895 

1 

V 

1.1283 75862 



e l = 

-7.8256 11091 

e 2 = 8.0068 96701 

e 3 * .38071 98964 

e 4 

.35851 70449 

e 5 * .079477 4487 



TABLE 17: 


FIFTH ORDER INTEGRATION CONSTANTS FOR 5 DERIVATIVE 
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7. FIFTH ORDER RUNGE-KUTTA TYPE INTEGRATOR 
WITH SIX DERIVATIVE EVALUATIONS 


As we saw in the previous section, we need more than five derivative 
evaluations to obtain a fifth order Runge-Kutta type integrator. In this 
section we attempt to obtain the fifth order integrator using six derivative 
evaluations, as shown below. 


*1 s 

ATf^ 

, t t + 

5-jAT) 






(76) 

II 

CSJ 

-*1 

ATfCx^ 

+ a ] k 1 

» t,| + ( 

$2^T ) 





(77) 

—3 = 

aT£(x. 

+ b^ 

+ bgkg 

, t^ + ( 

s 3 aT ) 




(78) 

II 

-T 

ATf(x. 

+ c l-l 

+ c z k 2 

+ c 3^3 

, ^ + , 

5 4 aT) 



(79) 

^5 = 

ATf(x. 

+ d ] k 1 

+ d 2 k 2 

+ d 3^3 

+ 

, t i + 

V T) 


(80) 

ii 

ATf(x. 

+ e l-l 

+ e 2-2 

+ c 3- k -3 

+ e ^4 

+ e 5^5 

* *1 + 

s 6 lT) 

(81) 

-i+1 

= *i + 

f l— 1 + 

f 2^2 + 

f 3^3 + 

f A + 

f 5^5 + 

f 6*6 

+ 0(1T 7 ) 

(82) 


The constraint equations that the integration constants must satisfy are 


6, » 0 

(83) 

a, - « 2 

(84) 

b l + b 2 ■ 5 3 

(85) 

C 1 + c 2 + c 3 = 6 4 

(86) 
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d l + 

d 2 + d 3 * d 4 * 4 5 



(87) 

e l + 

e 2 + e 3 + e 4 + e 5 3 S 6 



(88) 

V 

f 2 + f 3 + f 4 * f 5 + f 6 

= 1 


(89) 

f 2 6 2 

+ f 3*3 + f 4 5 4 + f 5 ! 5 + 

f 6 s 6 ’ 

= 1 

(90) 

f 2 S 2 

+ Vl + f 4*4 + Vb + 

f 5 S 6 ' 

= 1 

(91) 

Vf 

* f 3 s 3 + f 4 S 4 + f 5 5 5 + 

V* ' 

■1 

(92) 

f 2 S 2 

+ f 3 s 3 + f 4 s 4 + f 5 S 6 + 

f 6 6 6 ’ 

. 1 
-r 

(93) 


b 2 6 2 f 3 + ^ c 2 6 2 + c 3 6 3^ f 4 + ^ d 2 5 2 + d 3 6 3 + d 4 6 4^ f 5 
+ (e 2 6 2 + e 3 6 3 + e 4 6 4 + e 5 6 5^ f 6 = (T 

b 2 6 2 f 3 6 3 + ^ c 2 6 2 + c 3 6 3^4 6 4 + ^ d 2 5 2 + d 3 6 3 + d 4 6 4^ f 4 6 5 
+ (e 2 6 2 + e 3 6 3 + e 4 « 4 + e 5 6g)fg6g = g (9») 

b 2 6 2 f 3 6 3 + (c 2 S 2 + c 3 S 3 )f 4 S 4 + < d 2 S 2 + d 3 S 3 * d 4 S 4 )f 5 S 5 

+ (e 2 s 2 t e 3 s 3 + e 4 s 4 + e 5 s 5 )f 6 s| ■ jj (96) 

(b 2 6 2 ) 2f 3 + (C 2 S 2 * c.4 3 ) 2 f 4 ♦ (d 2 J 2 + d 3 o_ • d 4 S 4 ) 2 f 5 

+ (e 2 « 2 + e 3 « 3 + e 4 s 4 + e 5 « 5 ) 2 f 6 ■ jj (97) 

«>2«i f 3 * < c 2 { 2 + C 3 { 3 )f 4 + (d 2 5 2 * d 3 { 3 * d 4 S 4 )f 5 

+ (e 2 s| + e 3 6 3 + e 4 4 4 + e 5 «|)f 6 ■ {7 (98) 
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b 2‘i f 3 s 3 + <«2«2 • =3 s | )f 4 s 4 + {i A + d 3 5 3 + d 4 5 4> f 5 6 5 
* («2«1 1 e 3 S 3 * % S 4 * e 5 S 5 )f 8 4 6 * T5 <»> 

b 2 4 2 f 3 + ^ c 2 6 2 + C 3 S 3^4 + ' d 2 5 2 * d 3 4 3 * d 4 4 4 ,f 5 

+ (e 2 «| + « 3 «| * e 4 s| + e 5 s|)f 6 - jj (100) 

b 2 5 2^ c 3 f 4 6 4 + d 3 f 5 6 5 + e 3 f 6 5 6^ + ^ c 2 & 2 + c 3 6 3^ d 4 f 5 6 5 + e 4 f 6 5 6^ 

+ (d 2 5 2 + d 3 5 3 4- d 4 4 4 ).- 5 f 6 S 6 ■ (101) 

b 2 5 2 ( c 3 f 4 + d 3 f 5 + e 3 fg) + (c 2 6 2 + c^Md^fg + e^fg) 

+ < d 2«2 + d 3 s 3 + d 4 s 4 )fe 5 f 6 * ll (102) 

b 2 6 2^ c 3 f 4 + d 3 f 5 " r e 3 f 6^ 6 3 + ^ c 2 6 2 + c 3 6 3^ d 4 f 5 + e 4 f 6^ 6 4 
4 (d 2 6 2 4 d 3 J 3 4 d 4 5 4 )e 5 f 6 « 5 = (103) 

b 2*2^ c 3^4 + d 3 f 5 + e 3 f 6^ + ^ c 2 4 2 + c 3 4 3^ d 4 f 5 + e 4 f 6^ 

4 (d 2 «| 4 d 3 S‘ 4 d 4 «2)e 5 f 6 = ^y (104) 

t> 2 S 2 c 3^ d 4 f 5 * e 4 f 6^ b 2 i 2 d 3 e 5 f 6 + ^ c 2°2 + c 3 6 3^ d 4 e 5 f 6 ’ T’tf ^'° 5 ^ 

There are 23 constraint equations in 27 unknowns. This apparently allows 
us to supply 4 more constraint equations of our own choice. I picked the 
following 4 additional constraint equations. 
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f 2 6 2 + f 3 6 3 + f 4 6 4 + f 5«5 + f 6 fi 6 * l 


f/ 2 * f 3 6® + f 4 6$ + f 5 d® ♦ f g 6* = } 


f 2 6 2 + f 3 6 3 + f 4 5 4 + f 5 6 5 + f 6 6 6 * 8 


f/ 2 ♦ f,6® ♦ f 4 «® + f 5 «| ♦ f 6 «| * l 


(106) 

(107) 

(108) 
(109) 


These additional constraint equations will cause x = f(t) to be integrated 
to ninth order accuracy, giving us a 5/9 order integrator. 

We now have 27 equations In 27 unknowns. We tried diligently to find a 
solution to these equations. We used three different computer programs, each 
with a different method, but were only able to obtain an approximate solution 
of the 27 equations. Equations (83) through (89) were solved exactly. Equations 
(90) through (109) were solved to an average la accuracy of + 1.2*10"®. The 
integration constants that we obtained are shown in Table 18. We then made a 
brief attempt to find any solution of the original 23 equations in 27 unknowns. 

We were unsuccessful in that attempt, but the effort that we expended was not 
as great as it might have been. Thus we are reasonably sure that the 27 equa- 
tions in 27 unknowns has .<o exact solution. However, we can only say that we 
suspect that the 23 equations in 27 unknowns has no solution. 

One final comment concerning a sixth order Runge-Kutta type integrator. 

We have ascertained that the sixth order integration of x = f(x, t) requires 9 
derivative evaluations, giving 38 constraint equations in 45 unknowns, a truly 
formidable problem which again may have no exact solution. 
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« 4 * .71385 17583 

a 1 * .Out 5763 74117 

b 1 * -29.265 36037 

c 1 * 19.658 98673 

d-| = -18.777 16241 
d 4 = .30148 96548 

e 1 = 13.387 20508 
e 4 = .12617 48124 

f ] * -3.5372 73077 
f 4 = .17209 54083 


2 * .0015763 74117 
6 5 « .61937 61366 

b 2 * 29,570 50411 
c 2 = -19.883 49164 
d 2 * 19.272 14754 

e 2 = -13.325 28099 
e 5 * .44371 94067 

f 2 = 3.6754 69584 
f 5 = .18162 04139 


fi 3 * .30514 37420 
fi 6 = .93389 85460 


c 3 = .93835 66728 
d 3 = -.17709 86412 

e 3 = .30208 02382 

f 3 = .34296 07771 


f 6 = .16512 68939 


rABLE 18: 5/9 ORDER INTEGRATION CONSTANTS FOR 6 DERIVATIVE 

evaluations (AFpmtflAYE solut tot 
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8. A SELF-STARTING INTEGRATOR USING 3f/3x 


■ 

In conjunction with integrating x = f(x, t), the equations defining the 
matrix, af/ax, are frequently available. A self-starting integrator making 
use of this matrix is shown below. 


Given 


Let 


Then 


Where 



the proof is 


— i+.5 = — i+.5 + ^ 


Therefore 


3 a . t3 _. 


a(* 1+ 5 , t,+ .6aT) . i(x j+5 , t| + .5iT) + — cAT J • x 1+ 5 + o(ir) 


So 


il + l ■ *1 * ‘Tx, ♦ I it 2 x, - l 4T 2 x 1+ _ 5 + OUT 5 ) 
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But 


iV. 5 - i, ♦ ^ * 0(lT3) 


SO 


i i+ l = ^ + at ^ ^ + ^ 'x t + * 0(aT 3 ) 


t5> 


The integration coefficients in a self-starting integrator utilizing 
3f/ax, are unique. Thus there is no opportunity to adjust integrator constants 
to reduce truncation error, as there is with the Runge-Kutta type of integrators. 
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9. ESTIMATION OF TRUNCATION ERROR 


It Is frequently desirable to be able to estimate Integrator accuracy 
over the Integration Interval. This Information might then be used to control 
the Integration step size. A logical approach to this problem would to be 
make one extra derivative evaluation (generate one more JO* and then use this 
extra value of k, combined with the other k's, to estimate the truncation 
error. This would be equivalent to taking a, say, third order integrator, 
making one more derivative evaluation (evaluate k^), and obtaining a fourth 
order integration. This would then be coa^ared with the third order answer 
to determine the truncation error. Unfortunately it is not possible to do 
this. From Table 6 it is seen that a third order Integrator must have 


h - 8 3 V 6 2 
2 * « 2 2 - 3 <$ 2 

From Table 11 it is seen that a fourth order integrator must have 

w _ 6 3 *3 ' *2 

“2 • 26^ r~-T6 2 

It is clearly seen that the two values of b 2 can never equal one another. 

Thus a third order integrator can not be made into a fourth order integrator. 
In other words, you can't make a silk purse out of a sow's ear. 

In the following estimation of truncation error, we will assume that 
three integrations are made when integrating from t^ to t. +1 : one integration 
from t. to t. + the result being integrated from t. + 5 to t^; and one 
integration from t. to t^ +1 directly. 
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Integration from t^ to t^ + 5 , for nth order integrator, nay be repre- 
sented by 


^i+.5 ° -^-i’V + 


— i+.5 + 


AT 


kn+1 


n^. 1 

where a T = and where e(AT/2) is the truncation error associated with 

the particular integration scheme being used. 

Assuring e remains relatively constant, integration from t- + 5 to *i+l 
is given by 


or 


A 

A 


— i+1 



*i + 


.5 AT] 



1 


A 

A 


—i+1 



( 110 ) 


A single integration from to t,^, would yield 


-i+1 * -i+1 + ^ Ajn+1 


Subtracting equation (110) f^om (111) gives 


( 111 ) 



( 112 ) 
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As a first approximation ah/ax * I, so 



013 ) 


In situations where x = 



a better approximation is 


ah 

a7 



0 


IaT/2 

I 


In this case 



(114) 
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Knowing e, and assuming e remains relatively constant from one cycle to 
the next, a value of aT for the next cycle may be chosen to give arty desired 
integration accuracy. For example, suppose that x * £x" [ xj]^ and the error 
in ^ (the position error)is to be controlled. From equation (110), and as- 
suming ah/ax = I, the position error in the next cycle is seen to be 



Let the magnitude of the allowable position error vector be 


. = Units of position (feet, radians, etc.) 

°p Units of integration time (secs, days, etc.) 


Then 


6 AT 
P 


i+1 



Substituting equation (113) into the above, and solving for aT^ gives 



(115) 


The above equation has been used very successfully with a fourth order integra- 
tion of an earth to Mars trajectory. Here the position vector consisted of the 
positions of the earth, moon, sun, and Mars with respect to the spacecraft. 
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It was interesting to note that the motion of the earth-moon system was the 
dominant factor in controlling aT during most of the voyage. This is something 
that many people might overlook, yet It was readily detected using equation 
(115). 

When operating with very small values of 6 p (and small Integration step 
sizej), there may be a very large loss of significant figures In the quantity 
ft - 5. Generally, Integrator equations can be written In the form 


- 1+1 


+ 


Thus the loss of significant figures is dimished if we let x - ft ■ ax - aft. 

The reader may object that. In these cases, computer roundoff error Is the 
dominant error source - not truncation error. He Is correct. In these situa- 
tions 6p will not have absolute control over the integrator accuracy. Howev.r, 
the aT generated by the preceding equations will be an indicator of the "dynamic 
activity" of the system. 

We are now led to the interesting question of what is the best aT to use 
to obtain maximum integration accuracy. Theoretically, integration accuracy 
improves as aT -*• 0. However, as a practical matter, computer roundoff error 
will, at some point, cause the integration error to start increasing again as 
aT -*• 0. Suppose that the computer stores position accurate to (feet, 
radians, km, etc.). Again we will assume X is being used as the output of the 
integration. As before, the truncation error in position for one integration, for 
ft, is 



Now let the truncation error be the same size as the computer roundoff error. 
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Thus 




Using equation (113) and solving aT| + 1 gives 


AT 


1+1 




1 

-|n+T 


A 

, A A 

I*, - *, 


AT- 


A 

Had we been saving x instead of x as the output of the 
would have obtained the equation 


AT 


i+1 ‘ 2 


(2 n+1 -2k 


A 

, A A 

1-1 ' -1 


1 

n+T 


at. 


Note that a5-j - ax^ should be used in place of x-j - 


A 

A 

*1 


016) 


integrator, then we 


(117) 

in the above equations. 
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10. COMPUTER ALGORITHM FOR RUNGE-KUTTA TYPE INTEGRATORS 


The computing algorithm for Runge-Kutta type Integrators is particularly 
simple. For convenience of notation, let 


a l * A 11 

t>i , bg = A21 » Ag 2 

C 1 * c 2* c 3 = A 31 * A 32* A 33 


etc. 


Then, in engineering notation. Nth order integration of x * f(x) is accomplished 
by 


— i+1 * -i 
k_i * aT £(^+i ) 


-i+1 = -i + A ll— 1 
k 2 = ATf(x. +1 ) 

-i+1 = -i + A 21-l + A 22— 2 


h 9 

iii+1 = + A N1— 1 + A N2— 2 + * * • + A NN^N 
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In order to express the preceding equations In FORTRAN notation , let 
(for I « 1,2,...M and J « 1,2 N) 


X(I) - x 1+1 

XX(I) ■ 

^1 

F(I) » f 

DT - 

&T 

A(I,0) ■ A^j 

x(l.l) 

•*1 

K(1.2) * kg 

K(I.N) - 

h 

T * t 

TI-T, 

D ( 1 ) = 5i 

D(2) = 

6 2 ...D(N) 

* 6 N 


Then the FORTRAN equations for integrating x = f(x) are 



DO 1 

I = 1 ,M 

1 

xx(i) 

= x(i) 


DO 2 

J = 1 ,N 

comment 

EVALUATE F(l) = FUNCTION OF X(I) 


DO 2 

I = 1 ,M 


K(I,J) 

* DT*F(I) 


X(I) = 

XX(I) 


DO 2 

L * 1 ,J 

2 

X(I) * 

X(I) + A(J,L)*K(I,L) 


T * T 

► DT 


N* h Order Integrator for x * f (x) 
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The FORTRAN equations for Integrating x » f(x, t) are 



TI = T 



DO 1 

I = l.M 

1 

xx(i) 

= X(I) 


DO 2 

J = 1,N 


T = TI 

+ D(0)*DT 

comment 

EVALUATE F(I) = FUNCTION OF X(I),T 


DO 2 

I = 1 ,M 


K(I,J) 

= DT*F(I) 


X(I) = 

XX(I) 


DO 2 

L = 1 ,J 

2 

X(I) = 

X(I) + A ( J ,t)*K(I,L) 


TI = T 

I + DT 


N th Order Integrator for x - f(x, t) 


Though the above algorithms have the advantage of simplicity, the reader 
should be aware that if maximum speed of computing is desired, then the algo- 
rithms should be rewritten so that the doubly subscripted variables, A(I,J) 
and K(I,J), are replaced by singly subscripted variables. Singly subscripted 
variables require less tima to locate in core than do doubly subscripted 
variables. 
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11. SOME EMPERICAL EVALUATIONS 


The truncation e»*ror terms, In the previously discussed Integrators, are 
extremely complicate*!. The only practical way to investigate integrator ac- 
curacy is by emperical methods. One way to do this is to pick an x » f(x, t) 
equation whose solution is known and whose form is similar to the actual equa- 
tion of interest. We have done this here in this section for an orbit determina- 
tion problem. Another way is to use Eq. (116) (or similar equation) to obtain 
the most accurate solution of x = f(x^ t) that the computer can give. Other 
solutions would then be compared with this one. 

In order to provide an illustrative example of some of the integrators 
discussed in this report, a satellite orbit determination problem was run for 10 
orbits. The dynamical equation of motion, for a spherica 1 gravitational field, 
is given by the nonlinear, vector differential equation 


R 



R is the vector from the center of mass of the attracting body to the satellite, 
pis the gravitational parameter. The elements of the state vector were taken to 
be x, y, z, x, y, z. The state vector equation of motion, k - fjx, t), is then 


x 

y 

z 

X 

y 

z 


X 

y 

z 

-px/[x 2 + y 2 + 
-uy/[x 2 + y 2 + 
- u z/[x 2 + y 2 + 


z 2 ] 

z 2 J 

z 2 ] 


1.5 

1.5 

1.5 
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The initial conditions were precisely set to give a circular orbit, inclined at 
45°, and with a period of exact'y 6144 seconds. The accuracies of the various 
integrators with the above equation are shown in Tables 20 and 21. The term 
"2.5 order Taylor series" needs explaining. If x and x are elements of the state 
vector, then a second order Taylor series solution of the vector differential 
equation is 


x i+l * x i + x i A ^ + 

2 

x i+l “ *1 + + *j4T /2 

However, one would rarely use the equations in a computer program In this form. 
As long as x. is being computed, one would also use It to "improve" the estimate 
cf x. + .j. Thus we define a 2.5 order Taylor series solution to be 


x i+l = x i + x i AT + V t2/2 + V l3/6 


x i+l = x i + x i AT + 


Similarly a 3.5 order Taylor series solution goes out to fourth order in posi- 
tion and thiro order in velocity. Likewise for 4.5 and 5.5 solutions. 

The term "4/6 RK Set #r refers to the Runge-Kutta integration constants 
in Table 15. They provide fourth order integration of x = f(x) and 
sixth order for x = f(t). The set of integration constants "4th RK: 6 2 * 6 3 = * 5 ’ 

d 3 * .5" are fourth order Runge-Kutta integration constants which came from 
Table 12, where d 3 was optimized to give the best average position error over 10 
orbits of integration. 

The set of integration constants "4th RK: 6^ = .15, 6 3 * .192" deserves 

special mention. This set of fourth order Runge-Kutta constants was obtained 
from Table 11 by optimizing the two free parameters, and « 3 . The optimization 
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was not done with the orbit determination problem studied here. They were opti- 
mized with an entirely different orbit and with a different gravity field (equa- 
tions of motion).* Yet they worked extremely well with our current problem, which 
indicates that they probably would work very mil with any general orbit determina- 
tion problem. Their accuracy was even better than the 5th or 5.5 order Taylor 
series solutions. This outstanding set of integration constants Is shown in 
Table 19 below. 



TABLE 19: OPTIMUM FOURTH ORDER RUNGE-KUTTA INTEGRATION 

constants for OrbiT determination 


* Wm. M. Lear, "Direct Integration of Orbital Equations Using a Fourth Order 
Rur.ge-Kutta Integrator", TRW Technical Report 20029-601 3-T0-00, 31 August 
1972. 
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Table 20 shows the position error at the end of 10 orbits. Table 21 shows 
the average position error during the 10 orbits. In general the results of the 
two tables are in agreement as far as ranking the Integrators for accuracy. How- 
ever, there are a few cases of disagreement. For an Integration step size of 
256 seconds, the Runge-Kutta-Gtll integration constants had a final position 
error of 1274 meters, lower than with aT * 128 seconds which was 2193 meters. 
However, the average position error in Table 21 gives a truer picture of the 
overall accuracy. He see that for aT = 256 seconds, the Runge-Kutta-Gill con- 
stants had a rather large average error and thus Indicated that the absolute 
position error was passing through an abnormally low minimum at the end of 10 
orbits. In other words, the average position error Is a better indication of 
overall accuracy. He have not shown tables of velocity accuracy since, without 
exception, all the velocity errors are approximately 1/1000 of the position 
errors . 

Figures 1 and 2 show how the integrator error "grows" with time for 
several of the integrators. He see that the errors are not necessarily always 
increasing with time, but may exhibit periods where the error may decrease with 
time. Note particularly the excellent performance of the Table 19 integrator 
constants shown in Figure 2. 

Figure 3 shows the accuracy of various orders of Taylor series integrators 
versus integration step size. The important thing shown in this figure is the 
desirability of the fourth order integrator. Notice that for a constant inte- 
grator accuracy of 100 meters, the fourth order integration step size can be 13 
times larger than the third order integration step. Yet when we go from fourth 
order to fifth order, we can Increase our integration step by only a factor of 
1.3, hardly worth the effort of evaluating the extra derivative term. Also if 
one examines Tables 20 and 21, we also reach the conclusion that, at least for 
orbit integration, the fourth order integrators are perhaps the most desirable 
all purpose integrators. 
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12. INTEGRATION USING ACCELEROMETER DATA 


12.1 INTRODUCTION 

In aerospace work we frequently have the equation of motion 

R = a + a (118) 

— -g -s 


where 


R = 



* V 1, 11 


a^ is acceleration di*-e to gravity and is an analytic function of position and 
time (if higner order gravitational harmonics are used). The function is 
sensed acceleration, sensed by accelerometers, and is due to external forces 
other than gravity. However, the output of the accelerometers is not ^ itself, 
but the integral of a^. Thus what we have available for integration of the equa- 
tions of motion is 

v(t) =y‘ t a s (i)dT (119) 

For "nor.destruct" accelerometer data, t Q is the time at which the system was 
turned on. For "destruct" accelerometer data the integral is only over small 
time interval, say from t^ to We will assume nondestruct data here, the 

resulting equations can be easily modified for destruct data. 
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12.2 ONE COMMON SOLUTION OF EQ. (118) 

Solve (Integrate) the free-flight equation 


R = a^R, t) 


to obtain 


Rj^-j — f(Rj * Rj * t^ » AT) 


Rj + 1 * £,( Rj » Rjt ^ » ^T) 


Let 


**c 
P = —Js 

- 3R 


Note that P is a 3 by 3 matrix which is a function of £ and t. 
Then the solution of Eq. (118), 


R » ig + a^t) 


is given by 


Si*i 


t i+l 


f 11 

f(R f , Rj. t,, *T) *J SjfiHtj^j - t ) dx 


*! 




r^i+l 

S, +1 * a(«i. A,- V 1T > *1 


* 3T p iis > i‘ T3 + It [3P i^,i * * 


(120) 


( 121 ) 
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Let 


and note that 



a^-Odt * v^ +1 - exactly 
Thus our solution of the equations of motion Is given by 




( 122 ) 


( 123 ) 


( 124 ) 


( 125 ) 


Generally the last two t' ..is tne above equations are ignored. Note however, 
it is not particularly difficult to at least include the ^ «P^ (v. + j-v^ )/ aT 
term. Approximations for the integral, I_, are shown below. 
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AT/4 AT/4 aT/4 aT/4 

■< > — \ ■ ■ - i 

I * U C* 83 ^ + 32 li+i/4 + 1 ^-i+2/4 + 3 ^— 1+3/4 + 7 ^1+1 J 

f ! d Vi AT 7 



We have shown the higher order Integrators above, not because of their accuracy, 
but because the engineer doing the integration frequently has no control over 
when the accelerometers are read out. ine approximations of X shown aoove are 
much preferred to those shown below because of such things as engine cutoff in 
the middle of an Integration step. 
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AT 
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12.3 RUNGE-KUTTA INTEGRATION OF ACCELEROMETER DATA 

Runge-Kutta Integrators require the value of a^t) at various times within 
the integration Interval from t^ to t^ . The equations below supply this in- 
formation. Note that Runge-Kutta integration of accelerometer data is generally 
preferable to the previous method since the Runge-Kutta method easily accounts 
for the higher order tem.s, generally left out in the previous me "hod (Eqs. 

(124) and (125)). 




For 2nd Order R-K Integration 
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AT/4 aT/4 aT/4 aT/4 

| >- i 1 \ 




'i+1 


^*1 + «aT) = A^ + + A 2V i+<5 + A^.^g '■• AjV-^ 


. -5126 3 + 8646 ^ - 4166 + 48 

A 1 3aT 


. _ 7686 3 - 11526 2 + 4566 - 36 
A 2 = 3aT 


A 3 


-51 26 3 + 6726 2 - 2246 * 16 
3aT 


A 4 


1286 3 - 1446 2 +446-3 

_T 
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+ 6aT) = + A 2V i+>4 + AjV^g + A^v. + 8 + AgV i+1 


. _ 15 6256 4 - 35 0006 3 + 26 6256 2 - 77006 + 600 

A i zsn 


. _ -31 2506 4 + 65 0006 3 - 44 2506 2 + 10 7006 - 600 
A 2 24aT 


. _ 31 2506 4 - 60 0006 3 + 36 7506 2 - 78006 + 400 

*3 24aT 


. _ -15 6256 4 + 27 5006 3 - 15 3756 2 + 30506 - 150 
A 4 53aF 


. _ 31256 4 - 50006 3 + 26256 2 - 5006 + 24 
A 5 24aT 


A^ = -A,-A o -A 0 -A.-A c 
o 1 2 3 4 5 


da 5 

h - TOT® ( - 3750{5 * 937554 - 8500s3 

dt 


+ 33756* - 5486 + 24) 
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Note that any of the previous approximations of may be used with any lower 
or higher order integrator. The frequency with which the accelerometers are 
read out may determine the approximation that is used. The approximations of a^ 
shown above are much preferred to those shown below because of such things as 
engine cutoff in the middle of an integration step. 
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for 3rd Order R-K Integration 
(Accelerometer Sample Interval = <.T) 


AT AT aT aT 

I 1 y i 1 

t i-3 *i-2 t i -1 l i t i+l 

(t, + 6&T ) = A -3^i-3 + A _2— i-2 + A -l-i-l + A o^1 + A l-i+l 

_ 26 3 + 36 2 - 6 - 1 . -86 3 - 18A 2 +46+6 

A ~3 rST -2 s T2Z7 

. _ 126 3 + 366 2 +66-18 a _ 26 3 + 96^ +116+3 

A -l T2 aT A l V2 aT 

A o * * A -3" A -2' A -r A l 

H56 4 - S06 3 - 456 2 + 30o + 18) 

For 4th Order R-K Integration 
(Accelerometer Sample Interval = aT) 
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aT 


LT 


AT 


at 


t i-4 t i*3 t i-2 t i-l t i *i+l 

5s(tl + 6aT ^ = A -4^1-4 + A -3*i-3 + A -2^i-2 + A -l-i-l + A o^i + A l-i+l 


-56 4 - 206 3 - 15A 2 + 106+6 


256 4 + 1206 3 + 1056 2 - 606 - 40 


A -2 " 


-506 4 - 280 6 3 - 3306 2 + 1406 + 120 


506 4 + 3206 3 + 51 06 2 - 406 - 240 


5s 4 + 406 3 + 1056 2 + 1006 + 24 


A o = * A -4" A -3 _A -2' A -r A l 


4 V 5 , - , , 

8 - a ~ jj* 1 VW (- 6 « - 456^ - loos' 3 - 455* + 526 + 24) 
dt 3 
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APPENDIX A 


THE k EXPANSIONS 


Consider the vector differential equation, x * f(x). Let f 1 be the i th 
element of f. Let fl be the partial derivative of f’ with respect to the 
element of x. Let 


= fLfV 


= f’.rff* 


f jkx fJfkft 


fi - 


G ' * f ] fJ k/ f! 


* f ] f k fkft 


fjfkfifm 


ft. e f^f k f £f,P 

jk£ m 




f!.fj fk f £fm 

jk i m 


* f ]k fJ t f n fkf " 


f 3 f L fkf£f ” 


pi _ xifj rk~£jn 

v T j T kJtV T 


Q 1 « R 1 * fJriAM 

jktm j k £ m 


where the repeated indices indicate summation. 

Let the A vector be composed of the elements A*; let the B vector be com- 
posed of the elements B^, etc. Then it can be shown that 


A-l 



x = 8 
x » C + D 

T- E + 3F + G + H 

T« I_+6J + 4K+3L + 4M + N + 3P + a+ R 


1«i 

J 


ir f 1 A k 
jk e e 


ir f 1 A k e £ + 
5T f jkl e e E 


* + C 2 + C 3 

* d-j + d 2 + d 3 + 

* e l + e 2 * e 3 + e 4 + e 5 
ten t will get "bumped" by a-,&T, (b 


^+b 2 )AT > 


-2 








The fifth-order expansions of the k's are now shown below. 



-3 









k ^ « ATf (x + d 1 k 1 + dgkg + d^ + d^} * aTA + aT^jB + aT 3 6*C 
+ (d 2 * 2 + d 3 S 3 + d 4 6 4 )oj + aT 4 |^ 6^E + (d 2 6 2 + d 3 6 3 + d^JSgF 

+ -j(d 2 6 2 + d 3 « 3 + d 4 8 4 )G +■ ^2 5 2 4 3 + ^2 4 2 ^ ^3^3^ d 4^— ^ 

+ aT 5 |^ 6 &L + ^ d 2 6 2 + d 3 6 3 + d 4 6 4^ 5 5^ + ^ d 2 6 2 * d 3 6 3 + d 4 6 4^ 6 5^ 
+ + d 3 6 3 + d 4 5 4^ + ^2 fi 2 d 3 + ^ c 2^2 + ^3^3^ d 4^5~ 

+ ^(4^62 + ^3^3 + ^4^4)— * ^2 5 2 d 3 5 3 "*" ^ c 2 6 2 + c 3 6 3^ d 4 6 4^— 

■*■ ■j {b 2 62 d 3 * (^ 2^2 + c 3 fi 3^ d 4^ + ^2 fi 2 c 3 d 4^1 



kg » ATf(x + e 1 k 1 + e 2 k 2 + e^k 3 + + egkg) «* aTA + AT 2 6gB 


+ aT 6gC + ( e 2^2 + e 3 fi 3 + e 4 5 4 + e 5®5 


>a] + al4 [j 




* < e 2 ! 2 * e 3 S 3 + e 4 s 4 ♦ V 5 ) 5ft F + \ (e 2 «| + ^ * yj ♦ e g 6|)G 


* ^ 2 fi 2 e 3 * (c 2 6 2 * c 3 ® 3) e 4 * (^ 2 5 2 * ^363 + ^ 4 ^ 4 ) 6^}^ | 
+ A { &i * I < e 2 5 2 + e 3®3 + e 4 s 4 + e S S 5 )s |i 


+ \ <«2«| ♦ *3 S 3 * e 4 5 4 + *5 4 5> S <£ + 1 (e 2 s 2 * e 3 5 3 + e 4 { 4 + e S s S )2 t 


+ ■ 2® 3 ^ (c 2 6 2 ^ c 3 6 j ) ^ (d 2 6 2 + djjj + d^6^)e^}6gM 

+ i ^ e 2 6 2 + e 3 6 3 + e 4 6 4 + e 5 { 5^- + {b 2 6 2 e 3 6 3 + ^ c 2 6 2 + c 3 6 3^ e 4 6 4 

+ ( d 2 6 2 + d 3 6 3 + d 4 6 4 )e 5 6 5 } ^ + \ {b 2 6 2 e 3 + {c A + c 3 6 3 )e 4 

+ ( d 2<5 2 + d 3^3 + d 4 5 4^ e 5^3. * ^ b 2 S 2 C 3 e 4 + b 2^2^3^5 + ^ c 2^2 + ^3^3^ d 4 e 5^— 



APPENDIX B 


A FIFTH ORDER R-K INTEGRATOR WITH 
FIVE DERIVATIVE EVALUATIONS? 


Since one Infrequently hears about al' ^fcd fifth order Runge-Kutta inte- 
grators using five derivative evaluations, * appeared worthwhile to try and 
show that no such solution exists for the vector differential equation 
x * f (x) or x * f (x, t) . 

The Integrator constant constraint equations are as shown below. 


5j * 0 

(B-l) 

*1 " *2 

(B-2) 

b 1 h b 2 - 63 

(B-3) 

C 1 + c 2 + C 3 = 6 4 

(B-4) 

d 1 + d 2 + d 3 + d 4 - 6 5 

(B-5) 

®1 + e 2 ♦ e 3 + e 4 + 65 - 1 

(B— 6 ) 

e 2 6 2 + e 3 6 3 + e 4 6 4 + e 5 6 5 " 

(B-7) 

e 2 5 2 + e 3 { 3 + e 4 6 4 + e 5 S 5 * 1/3 

(B- 8 ) 

e 2 { 2 + e 3 6 3 + *4 6 4 + e 5 fi 5 " ^ 

(B-9) 

®2 { 2 + e 3 6 3 + e 4 6 4 + e 5 6 5 “ 

(B-l 0) 


B-l 



b 2 d 2 e 3 + ^ c 2 5 2 * c 3 fi 3^ e 4 * ^ d 2^2 * ^3^3 + d 4^4^®5 “ 1/6 (B-ll) 

^2 fi 2 e 3 S 3 + ^ c 2 fi 2 + c 3 5 3^ e 4 5 4 + ^2^2 + ^3^3 + d 4^4^ e 5^5 * 1/® (B— 12) 

* > 2 <S 2 e 3 6 3 + ^ c 2 6 2 ^ c 3 5 3^ e 4 6 4 ^ d 2^2 ^3^3 + ^4^4^5^5 s 1/10 ( B — 13) 

(b 2 6 2 ) 2 e 3 + (c 2 6 2 + c 3 6 3 ) 2 e 4 + (d 2 6 2 + d 3 fi 3 + d 4 6 4 ) 2 e 5 * 1/20 (B-14) 

b 2 6^e 3 ♦ (c 2 6^ + c 3 6 3 )e 4 + (d 2 « 2 + d 3 6 2 + d 4 6 2 )e 5 = 1/12 (B-15) 

b 2 6 2 e 3 6 3 + (c 2 6 2 + c 3 6 3 )e 4 6 4 + (d 2 6 2 + d 3 6 3 + d 4 6 4 )e 5 6 b 1 ' (B-16) 

b 2 6 2 e 3 + (c 2 6 2 + c 3 6 3 )e 4 + (d 2 6 2 + d 3 6 3 + d 4 6 4 )e 5 = lA (B-17) 

b 2^2 C 3 e 4 5 4 + - b 2^2 d 3 + ^ c 2^2 + c 3^3^ d 4^ e 5^5 * ^/ BB (B-18) 

b 2 fi 2 C 3 e 4 + ^- b 2 <S 2 d 3 + ^2^2 + ^3^3^ d 4^®5 = 1 /^4 (B-19) 

b 2 6 2 c 3 e 4 5 3 + ^ b 2 6 2 d 3 S 3 + ^2^2 + ^3^3^ d 4^4^ e 5 = 1/40 (B-20) 

b 2 5 2 c 3 e 4 ^ b 2 S 2 d 3 ^ 2^2 * ^ 3 ^ 3 ^ d 4^®5 ~ 1/00 ( B — 21 ) 

b 2 6 2 c, 3 d 4 * 5 « 1/120 (B-22) 


22 equations In 20 unknowns. Solving Eqs. (B-7) through (B-10) gives 


B-2 



(« 4 + « 5 )(206 3 - 15) ♦ 5 4 6 5 (20 - 30fi 3 ) * 12 - 15« 3 

®2 M $ Ofi 2 (« 2 - « 3 )( fi 2 - « 4 )(« 2 - fig ) 


(B-23) 


(«4 + fig) (206 2 - 15) + fi 4 fig(23 - 30« 2 ) + 12 - 154 2 
606 3 i6 3 - 6 2 )(fi 3 ‘ 6 4)( fi 3 ~ *5) 


(B-24) 


e 4 


(6 2 + fi 3 )(206 5 - 15) + 6 2 6 3 (20 - 306 5 ) + 12 - 156 g 

6064(64 - fi 5 )(fi 4 - « 3 )(6 4 - 6 2 ) 


(B-25) 


e 5 = 


(6 2 + 6 3 )(206 4 - 15) + 6 2 5 3 (20 - 306 4 ) + 12 - 156 4 

= 606g(6g - 6 4 H«5 ~ « 3 )(fig - fi 2 ) 


(B-26) 


Equations (B-ll) and (B-12) give 


b 2 6 2 e 3^ 6 3 " + ^ c 2 6 2 + c 3 fi 3^ e 4^4 ~ = ~W 


3-46, 


(B-27) 


Equations (B-12) and (B-13) give 


4-56, 


b 2 6 2 e 3 6 3^ 6 3 ‘ 6 5^ + ^ c 2 d 2 + c 3 6 3^ e 4 6 4^ 6 4 “ 6 5^ = ~W 


(B-28) 


The above two equations yield 


3-3 



(B-29) 


-15(« 4 ♦ 6 g ) ♦ 20 6 4 6 5 + 12 

b 2 4 2 * --- s s r ' 


Substituting Eq. (B-24), for e^. Into Eq. (B-29) gives 


6 3 (« 3 - « 2 )[-15(6 4 + 6 g ) + 20 « 4 « 5 + 12] 

b 2 fi 2 * m 2 -■ so) c& 4 (wwftjrsjB"”* n - yk 2 


(B-30) 


Now Eqs. (B-19) and (B-21) yield 


. d _ _ 1 2 ” 56 2 
c 3 a 4 e 5 * 63 - « 2 


(B-31 ) 


But from Eq. (B-22) 


c 3 d 4 e 5 = m 


(B-32) 


Setting Eq. (B-31) equal to Eq. (B-32) gives 


^ 43(63 - 6 2 ) 

b 2 6 2 ' 2 - 56 2 


(B-33) 


Setting Eq. (B-30) equal to Eq. (B-33) yields 


6p[(7 - 8<$ 4 )6g + 7$ 4 - 6] = 0 


B \ 



But fron Eq. (B-22), 6 2 t 0, therefore 

7« 4 - 6 
s 5' 55^7 

Now Eqs. (B-19) and (B-20) give 

] 3 - 5i^ 

c 3 e 4 * d 3 e 6 * 1% b 2 «j « 3 - 8, 

Substituting Eq. (B-33) Into the above yields 

(2 - 56 2 )(3 - 56 4 ) 

C 3 e 4 + d 3 e 5 * 1*0 6"t« 3 - 6 2 )(« 3 "- « 4 ) 
Equations (B-17) and (B-15) give 

3 - ! 

c 3 6 3^ 6 3 ~ 6 2^ e 4 + d 3 6 3^ 6 3 " 6 2^ e 5 + d 4 6 4^ 6 4 ' 6 2^ e 5 * “ S5” 
Equations (B-17) and (B-ll) give 

3 _ 

0363(63 - 6 2 )e 4 + ^363(63 - 6^)e 5 + d 4 6 4 (6| - 6^ * — 


(B-34) 


(B-35) 


(B-36) 


(B-37) 


(B-38) 


B-5 



Combining the above two equations yields 


c 3 e 4 + d 3 e 5 “ 


- 10« 4 - 106 2 + 206 2 6 4 + 6 


> 3'°3 " °2 ;v 3 


Setting this equation equal to Eq. (B-36) results in 


(B— 39) 


« 2 (« 4 - 1 ) * 0 

But from Eq. (o-22) we see that « 2 t 0. Therefore 

5 4 = 1 (B— 40) 

Substituting this into Eq. (B-34) yields 

6 5 = 1 (B-41 ) 

Now comparing Eq. (B-18) with (B-19), we see that i 4 Mj ! 1 can not be a 
solution. Thus we have a contradiction, and there appears to be no valid 
solution. Also in going back to investigate the singular points of the various 
equations, I still find contradictions. And subjecting the above equations to 
two different least-squares, computer solutions, I find no exact solutions. 

Thus it appears that one more derivative evaluation is needed for a fifth order 
Runge-Kutta type of integrator for the vector differential equation, x * f(x, t). 


8-6 



